%% Navier-Stokes equation solver
%
%
%

%% Problem settings

global NIt
global NJt
global NIe
global NJe

% Domain dimensions

Lxt  = 12;
Lxe  = 4;       % Longitude before expansion
Lyt  = 2;
Lye  = 1;

% Discretization, pressure finite volumes

% I denotes x direction
% J denotes y direction

elemsPerLenghtUnit = 6;

NIt = Lxt*elemsPerLenghtUnit;
NIe = Lxe*elemsPerLenghtUnit;
NJt = Lyt*elemsPerLenghtUnit;
NJe = Lye*elemsPerLenghtUnit;

% Nodal spacing

h   = Lxt/NIt;

% Constitutive parameters

nu  = 0.1;
rho = 1;

% Boundary conditions


% Time domain and discretization
dt = 0.05;
T0 = 0;
T1 = 40.;
nSkip = 1;             % skips saving solution
nTimeSteps = (T1-T0)/dt;
nSaveSteps =  (T1-T0)/(dt*nSkip);
nSave      = 0;

dtLim1 = 0.25 * h^2 / nu
dtLim2 = 2*nu/1

if (dt>dtLim1 || dt>dtLim2 )
fprintf ('stability warning')
end



%% Pre-processing
%  Memory allocation

dimP = NIt*NJt-NIe*(NJt-NJe);
dimU = dimP + NJt;
dimV = dimP + NIt;

% Problem variables
P = zeros(dimP,1);
U = zeros(dimU,1);
V = zeros(dimV,1);

% Auxiliars: midstep variables & Pk for iterative pressure scheme computation
Pk = zeros(dimP,1);
Us = zeros(dimU,1);
Vs = zeros(dimV,1);

% Saving solution over time
PT = zeros(dimP,nSaveSteps);
UT = zeros(dimU,nSaveSteps);
VT = zeros(dimV,nSaveSteps);

% Auxiliar linear matrix & load vector for pressure implicit equation
M = zeros (dimP,dimP);
l = zeros (dimP,1);


% Dimensions

NIut = NIt+1;
NJut = NJt;
NIue = NIe;
NJue = NJe;

NIvt = NIt;
NJvt = NJt+1;
NIve = NIe;
NJve = NJe+1;

% Mapping spatial points

Xp = PspatialMap( h , Lxt , Lxe , Lye , Lyt , NIt , NJt , NJe , NIe , dimP);
Xu = UspatialMap( h , Lxt , Lxe , Lye , Lyt , NIut , NJut , NJue , NIue , dimU);
Xv = VspatialMap( h , Lxt , Lxe , Lye , Lyt , NIvt , NJvt , NJve , NIve , dimV);

assemblePressureMatrix;



%% Processing - Predictor step
%
%
%
%

% Temporal loop

time = T0;
nstep = 0;
while (time<T1)
    
    
    
    time = time + dt;
    nstep = nstep+1;
    
    fprintf('Solving timestep %d of %d \n',nstep,nTimeSteps)
    
    % Block 1
    
    % Velocity in X
    for i=1:NIut
        for j=1:NJue
            
            U_cen   = U ( Umap (i  ,j  ) );  % center of finite volume
            
            % Boundary conditions on inlet & outlet
            if (i==1) Us ( Umap (i,j) ) = getParabola ( Xp ( Umap (i,j),2 ) ) ; continue ; end    % Boundary condition inlet
            if (i==NIut) Us ( Umap (i,j) ) = Us ( Umap (i-1,j) ) ; continue ; end    % Boundary condition oulet - null normal gradient
            
            % Right and left volumes
            U_right = U ( Umap (i+1,j  ) );
            U_left  = U ( Umap (i-1,j  ) );
            
            % Top and bot volumes, considering ghost nodes
            if ( j==NJue && i<=NIue ) U_top   = - U_cen; else U_top   = U ( Umap (i  ,j+1) ); end
            if ( j==1 ) U_bot = - U_cen ; else U_bot   = U ( Umap (i  ,j-1) ); end
            
            % Velocity in Y direction on vertices of the finite volume
            V_tr    = V ( Vmap (i  ,j+1) );  % top-right
            V_tl    = V ( Vmap (i-1,j+1) );  % top-left
            V_br    = V ( Vmap (i  ,j  ) );  % bot-right
            V_bl    = V ( Vmap (i-1,j  ) );  % bot-left
            
            Ax  = + ( U_cen + U_right )^2 ...               % Right boundary
                + ( V_tr + V_tl ) * ( U_cen + U_top ) ...   % Top boundary
                - ( U_cen + U_left  )^2 ...                 % Left boundary
                - ( V_br + V_bl ) * ( U_cen + U_bot );      % Bottom boundary
            
            Ax  = Ax/(4*h);
            
            Dx  = U_left + U_right + U_top + U_bot - 4 * U_cen;
            
            Dx  = Dx/h^2;
            
            Us ( Umap (i,j) ) = U_cen + dt* (-Ax+nu*Dx);
            
        end
    end
    
    % Velocity in Y
    for i=1:NIvt
        for j=1:NJve
            
            V_cen  = V ( Vmap (i  ,j  ) );  % center of finite volume
            
            % Boundary conditions on top and bot boundaries
            if (j==1) Vs (Vmap (i,j)) = 0; continue ; end
            if (j==NJve && i<=NIve ) Vs (Vmap (i,j)) = 0; continue ; end
            
            % Right and left boundaries, considering ghost nodes
            if (i==1) V_left = -V_cen; else V_left = V ( Vmap (i-1,j ) ); end
            if (i==NIvt) V_right = -V_cen; else V_right = V ( Vmap (i+1,j ) ); end
            
            V_bot = V ( Vmap (i,j-1));
            V_top = V ( Vmap (i,j+1));
            
            U_tr    = U ( Umap ( i+1,j  ) );
            U_tl    = U ( Umap ( i  ,j  ) );
            U_br    = U ( Umap ( i+1,j-1) );
            U_bl    = U ( Umap ( i  ,j-1) );
            
            Ay  = + ( U_tr + U_br ) * ( V_cen + V_right ) ...
                + ( V_cen + V_top )^2 ...
                - ( U_tl + U_bl ) * ( V_cen + V_left ) ...
                - ( V_cen + V_bot )^2;
            
            Ay  = Ay/(4*h);
            
            Dy  = V_left + V_right + V_top + V_bot - 4* V_cen;
            
            Dy  = Dy/h^2;
            
            Vs ( Vmap (i,j) ) = V_cen + dt* (-Ay+nu*Dy);
            
        end
    end
    
    
    % Block 2
    
    % Velocity in X
    for i=NIue+1:NIut
        for j=(NJue+1):NJut
            
            U_cen   = U ( Umap (i  ,j  ) );  % center of finite volume
            
            % Boundary conditions on expansion & outlet
            if (i==NIue+1) Us ( Umap (i,j) ) = 0 ; continue ; end    % Boundary condition expansion
            if (i==NIut) Us ( Umap (i,j) ) = Us ( Umap (i-1,j) ) ; continue ; end    % Boundary condition oulet - null normal gradient
            
            % Right and left volumes
            U_right = U ( Umap (i+1,j  ) );
            U_left  = U ( Umap (i-1,j  ) );
            
            % Top and bot volumes, considering ghost nodes
            if ( j==NJut ) U_top   = - U_cen; else U_top   = U ( Umap (i  ,j+1) ); end
            
            % Velocity in Y direction on vertices of the finite volume
            V_tr    = V ( Vmap (i  ,j+1) );  % top-right
            V_tl    = V ( Vmap (i-1,j+1) );  % top-left
            V_br    = V ( Vmap (i  ,j  ) );  % bot-right
            V_bl    = V ( Vmap (i-1,j  ) );  % bot-left
            
            Ax  = + ( U_cen + U_right )^2 ...
                + ( V_tr + V_tl ) * ( U_cen + U_top ) ...
                - ( U_cen + U_left  )^2 ...
                - ( V_br + V_bl ) * ( U_cen + U_bot );
            
            Ax  = Ax/(4*h);
            
            Dx  = U_left + U_right + U_top + U_bot - 4 * U_cen;
            
            Dx  = Dx/h^2;
            
            Us ( Umap (i,j) ) = U_cen + dt* (-Ax+nu*Dx);
            
        end
    end
    
    % Velocity in Y
    
    for i=NIve+1:NIvt
        for j=NJve+1:NJvt
            
            V_cen  = V ( Vmap (i  ,j  ) );  % center of finite volume
            
            % Boundary conditions on top and bot boundaries
            if (j==NJvt ) Vs (Vmap (i,j)) = 0; continue ; end
            
            % Right and left boundaries, considering ghost nodes
            if (i==NIve+1) V_left  = -V_cen; else V_left  = V ( Vmap (i-1,j ) ); end
            if (i==NIvt)   V_right = -V_cen; else V_right = V ( Vmap (i+1,j ) ); end
            
            V_bot = V ( Vmap (i,j-1));
            V_top = V ( Vmap (i,j+1));
            
            U_tr    = U ( Umap ( i+1,j  ) );
            U_tl    = U ( Umap ( i  ,j  ) );
            U_br    = U ( Umap ( i+1,j-1) );
            U_bl    = U ( Umap ( i  ,j-1) );
            
            Ay  = + ( U_tr + U_br ) * ( V_cen + V_right ) ...
                + ( V_cen + V_top )^2 ...
                - ( U_tl + U_bl ) * ( V_cen + V_left ) ...
                - ( V_cen + V_bot )^2;
            
            Ay  = Ay/(4*h);
            
            Dy  = V_left + V_right + V_top + V_bot - 4* V_cen;
            
            Dy  = Dy/h^2;
            
            Vs ( Vmap (i,j) ) = V_cen + dt* (-Ay+nu*Dy);
            
        end
    end
    
    
    
    
    
    %% Processing - pressure step
    %
    
    c2 = rho/(dt*h);
    
    for i=1:NIt
        for j=1:NJt
            
            if (i<=NIe && j>NJe) continue; end
            
            U_right = Us ( Umap (i+1,j  ) );
            U_left  = Us ( Umap (i  ,j  ) );
            V_top   = Vs ( Vmap (i  ,j+1) );
            V_bot   = Vs ( Vmap (i  ,j  ) );
            
            if (i~=NIt)
                l ( Pmap (i,j) ) = U_right - U_left + V_top - V_left;
            else
                l ( Pmap (i,j) ) = 0 ;
            end
            
        end
    end
    
    l = l/c2;
    %     if (nstep==1)
    %         M(1,:) = 0;
    %         M(1,1) = 1;
    %         l (1) = 0;
    %     else if (nstep==2)
    %             assemblePressureMatrix;
    %         end
    %     end
    
    P = M\l;
    
    
    
    %% Processing - projection step
    
    c1 = (dt/(rho*h));
    
    % Velocity X
    
    for i=1:NIut
        for j=1:NJut
            if (i<=NIue && j>NJue) continue; end
            
            % Top boundaries
            if ( j==NJut ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            if ( (j==NJue) && (i<=NIue) ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            
            % Bot boundary
            if ( j==1 ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            
            % Right boundary
            if ( i==NIut ) U ( Umap (i,j) ) = U ( Umap (i-1,j) ) ; continue ; end
            
            % Left boundary
            if ( i==1 ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            if ( i==NIue+1 && j>NJue ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            
            
            
            U ( Umap (i,j) ) = Us ( Umap (i,j) ) - c1 * ( P ( Pmap(i,j) ) - P ( Pmap(i-1,j) ) ) ;
        end
    end
    
    % Velocity Y
    
    for i=1:NIvt
        for j=1:NJvt
            if (i<=NIue && j>NJue) continue; end
            
            % Top boundaries
            if ( j==NJvt ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            if ( (j==NJve) && (i<=NIve) ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            % Bot boundary
            if ( j==1 ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            % Right boundary
            if ( i==NIvt ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            % Left boundary
            if ( i==1 ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) - c1 * ( P ( Pmap(i,j) ) - P ( Pmap(i,j-1) ) );
        end
    end
    
    % Saving results
    
    if ( mod(nstep-1,nSkip) == 0 )
        nSave = nSave+1;
        PT (:,nSave)= P;
        UT (:,nSave)= U;
        VT (:,nSave)= V;
        
    end
    
end


%% Post-Process



% Plot pressure
figure()

scatter3 (Xp(:,1),Xp(:,2),PT(:,10),'r','.')
title('pressure');

% Plot velocity x
figure()
%for i=1:1000
scatter3 (Xu(:,1),Xu(:,2),UT(:,10),'k','.')
title('velocity u');

%end
% Plot velocity y
% figure()
%
% scatter3 (Xv(:,1),Xv(:,2),VT(:,750),'g','.')
% title('velocity v');






























